SKI WASHING

1 Introduction

Located in the Swiss Alps, the ski resorts Engelberg-Titlis and Rigi-Kaltbad, along with the communities of Engelberg and Küssnacht, are facing significant challenges due to climate change. Especially winter tourism is crucial for these areas, contributing 1% of Switzerland’s GDP and 10-20% of economic output in certain regions. The increasingly unpredictable reliability of snowfall, undermines the examination of snowmaking technologies to sustain ski resort operations and regional economies.

The analysis is guided by the principles of Sustainable Development Goals 7 Affordable and Clean Energy, 8 Decent Work and Economic Growth, and 9 Industry, Innovation, and Infrastructure, focusing on energy efficiency, economic sustainability, and innovation in addressing climate change impacts on winter tourism.

Considering th various stakholders involved, including Titlis and Rigi Bergbahnen, the International Ski Federation (FIS), TechnoAlpin, public transportation, and local culture, the analysis adopts a comprehensive approach. This investigation encompasses meteorological and snow cover data from 1990 to 2024, alongside guest occupancy rates from 2013 onwards. By applying this method the future viability of snowmaking as a method to mitigate climate change effects on ski resorts shall be assessed.

The objective is to evaluate sustainability of snowmaking technology in maintaining ski resort operations against changing climate conditions, offering insights into the resilience and adaptability of the winter tourism in Switzerland.

1.1 Research Question

In a scientific pursuit to inform strategic foresight for ski resorts, our core research question is articulated as follows: “How might we predict the duration of effective snow-making in coming years using weather trends, and understand its consequences for winter sports destinations?”

1.2 Method

Our investigation employs time series analysis using R to dissect decades of weather and snow data, alongside tourism occupancy rates. Key R packages aid in data cleansing and trend decomposition, revealing patterns and cyclicality. We perform correlation analysis to uncover relationships between snowmaking, weather, and tourism. The predictive phase uses the ARIMA model to forecast future snowmaking needs and the sustainability of ski seasons, providing data-driven insights for strategic planning and long-term resort viability.

1.3 Assumptions and Limitations

This study is conducted with the acknowledgment of several constraints. Notably, gaps exist within our datasets, with certain years missing valuable data points, and location-specific weather information proving elusive. Additionally, hotel occupancy rates, a key indicator of economic vitality, are available only from 2013 forward. These limitations necessitate assumptions to bridge data voids and infer trends. While such estimations enable a continuous analytical narrative, they also introduce a degree of uncertainty. As such, the findings and predictions presented must be interpreted with an understanding of these inherent constraints.

2 EDA

2.1 Temperature from homogenous data source

homogenous data: Adjusting historical measurements to current standards

# Laden der installierten Pakete
library(ggplot2)
library(dplyr)
library(stats)
library(zoo)

# Pfad zum Unterordner
data_folder <- "Daten/"

# Funktion zum Extrahieren des Standortnamens aus dem Dateinamen
get_location <- function(file_name) {
  if (grepl("Engelberg", file_name)) {
    return("Engelberg")
  } else if (grepl("Luzern", file_name)) {
    return("Luzern")
  } else {
    return("Unknown")
  }
}

# Laden der Datensätze für Engelberg und Luzern
Engelberg <- read.csv(paste0(data_folder, "Messreihe_Engelberg.csv"))
Luzern <- read.csv(paste0(data_folder, "Messreihe_Luzern.csv"))

# Extrahieren des Standortnamens und Hinzufügen zur Datensätze
Engelberg$Location <- get_location("Messreihe_Engelberg.csv")
Luzern$Location <- get_location("Messreihe_Luzern.csv")

# Zusammenführen der Daten
combined_data <- rbind(Engelberg, Luzern)

# Umwandlung der Spalten 'Year' und 'Month' in ein Datumsformat
combined_data$Date <- as.Date(paste(combined_data$Year, combined_data$Month, "01", sep = "-"))

# Gruppieren nach Jahr und Standort und Berechnen des Durchschnitts pro Jahr
summary_data <- combined_data %>%
  group_by(Year, Location) %>%
  summarise(Avg_Temperature = mean(Temperature),
            Avg_Precipitation = mean(Precipitation))

# Plotten der Temperaturzeitreihe für Engelberg und Luzern
avg_temp_plot <- ggplot(summary_data, aes(x = Year, y = Avg_Temperature, color = Location)) +
  geom_line() +
  labs(x = "Year", y = "Average Temperature (°C)", color = "Location") +
  ggtitle("Average Temperature Time Series for Engelberg and Luzern") +
  theme_minimal()

# Plotten der Niederschlagszeitreihe für Engelberg und Luzern
avg_preci_plot <- ggplot(summary_data, 
                         aes(x = Year, y = Avg_Precipitation, color = Location)) +
  geom_line() +
  labs(x = "Year", y = "Average Precipitation (mm)", color = "Location") +
  ggtitle("Average Precipitation Time Series for Engelberg and Luzern") +
  theme_minimal()

# Gruppieren nach 5-Jahres-Schritten und Berechnen des Durchschnitts pro Zeitraum
summary_data_5yr <- combined_data %>%
  mutate(Year_Group = as.integer((Year - 1) %/% 5) * 5 + 1) %>%
  group_by(Year_Group, Location) %>%
  summarise(Avg_Temperature = mean(Temperature),
            Avg_Precipitation = mean(Precipitation))

# Plotten der Temperaturzeitreihe für Engelberg und Luzern (5-Jahres-Periode)
avg_5temp_plot <- ggplot(summary_data_5yr, aes(x = Year_Group, y = Avg_Temperature, color = Location)) +
  geom_line() +
  labs(x = "5-Year Period", y = "Average Temperature (°C)", color = "Location") +
  ggtitle("Average Temperature Time Series for Engelberg and Luzern (5-Year Period)") +
  theme_minimal()

# Plotten der Niederschlagszeitreihe für Engelberg und Luzern (5-Jahres-Periode)
avg_5preci_plot <- ggplot(summary_data_5yr, 
                          aes(x = Year_Group, y = Avg_Precipitation, color = Location)) +
  geom_line() +
  labs(x = "5-Year Period", y = "Average Precipitation (mm)", color = "Location") +
  ggtitle("Average Precipitation Time Series for Engelberg and Luzern (5-Year Period)") +
  theme_minimal()

# Konvertieren der aggregierten Daten in eine Zeitreihe
temp_ts <- ts(summary_data$Avg_Temperature, start = c(1864, 1), frequency = 12)
precip_ts <- ts(summary_data$Avg_Precipitation, start = c(1864, 1), frequency = 12)

# Fehlende Werte in der Zeitreihe ersetzen (hier mit 0)
temp_ts_filled <- na.fill(temp_ts, 0)
precip_ts_filled <- na.fill(precip_ts, 0)

# Anwenden der decompose-Funktion
temp_decomposed <- decompose(temp_ts_filled)
precip_decomposed <- decompose(precip_ts_filled)

2.1.1 Visualizing

We will be using the homogenous data from Luzern to compliment the Seebodenalp (Rigi) temperature data.

avg_temp_plot

The plot displays a time series of average temperatures for Engelberg and Luzern, with the x-axis representing the years from approximately 1900 to just beyond the year 2000 and the y-axis indicating temperature in degrees Celsius. Engelberg’s temperature is depicted with a red line, while Luzern’s temperature is shown with a blue line. Both locations show fluctuations that suggest seasonality and an overall upward trend in temperatures over the century, indicating a rise in average temperatures.

avg_5temp_plot

The plot illustrates a time series of average temperatures over a 5-year period for Engelberg (red line) and Luzern (blue line), showing a clear upward trend in temperatures for both locations from 1900 to the 2000s.

avg_preci_plot

The graph shows the average precipitation time series for Engelberg and Luzern, with both locations showing significant variability in precipitation over the years from 1900 to beyond 2000. The Engelberg data is depicted with a red line and exhibits more pronounced spikes, whereas the Luzern data, represented by a blue line, shows somewhat less volatility in average precipitation levels.

avg_5preci_plot

# Plotten der saisonalen, trendigen und saisonal bereinigten Komponenten der Temperaturzeitreihe
plot(temp_decomposed)

The plot presents a temperature decomposition of an additive time series from 1865 to 1890 into three components: observed, trend, and seasonal. The observed component shows the actual time series data with fluctuations, while the trend component indicates a relatively stable or slightly increasing trend over the 25-year period. The seasonal component displays consistent, periodic fluctuations suggesting a regular pattern that repeats annually, while the random component, which represents the residual noise, shows irregularities that cannot be attributed to the trend or seasonal effects.

# Plotten der saisonalen, trendigen und saisonal bereinigten Komponenten der Niederschlagszeitreihe
plot(precip_decomposed)

The chart presents a preciption decomposition of an additive time series from 1865 to 1890 into three components: observed, trend, and seasonal. The observed component shows the actual time series data with fluctuations, while the trend component indicates a relatively stable or slightly increasing trend over the 25-year period. The seasonal component displays consistent, periodic fluctuations suggesting a regular pattern that repeats annually, while the random component, which represents the residual noise, shows irregularities that cannot be attributed to the trend or seasonal effects.

2.2 Temperature from meteo swiss

library(tidyverse)

#### Daten einlesen ####

# Engelberg
engelberg <- read.csv("Daten/Engelberg_data.csv", sep = ';', na.strings = c('-',''))

engelberg$time <- as.Date(as.character(engelberg$time), format = "%Y%m%d")
engelberg_clean <- engelberg %>% select('stn',
                                        'time',
                                        'tre200nn')

# Seeboden
seeboden <- read.csv("Daten/Seebodenalp_data.csv", sep = ';', na.strings = c('-',''))
# Convert integer dates to Date objects
seeboden$time <- as.Date(as.character(seeboden$time), format = "%Y%m%d")
seeboden_clean <- seeboden %>% select('stn',
                                      'time',
                                      'tre200nn')

Engelberg: inspection of the needed variables

str(engelberg_clean)

## 'data.frame':    15737 obs. of  3 variables:
##  $ stn     : chr  "ENG" "ENG" "ENG" "ENG" ...
##  $ time    : Date, format: "1980-01-01" "1980-01-02" ...
##  $ tre200nn: num  -5.6 -6.8 -6.5 -9.8 1 -3.4 -1.4 -8.4 -12.6 -4.5 ...

summary(engelberg_clean)

##      stn                 time               tre200nn      
##  Length:15737       Min.   :1980-01-01   Min.   :-25.100  
##  Class :character   1st Qu.:1990-10-09   1st Qu.: -1.500  
##  Mode  :character   Median :2001-07-17   Median :  3.400  
##                     Mean   :2001-07-17   Mean   :  3.079  
##                     3rd Qu.:2012-04-24   3rd Qu.:  8.600  
##                     Max.   :2023-01-31   Max.   : 21.400  
##                                          NA's   :654

Seeboden: inspection of the needed variables

# NA anschauen
str(seeboden_clean)

## 'data.frame':    11889 obs. of  3 variables:
##  $ stn     : chr  "NABRIG" "NABRIG" "NABRIG" "NABRIG" ...
##  $ time    : Date, format: "1991-06-14" "1991-06-15" ...
##  $ tre200nn: num  9.4 11.3 12.7 6.2 NA 5 6.2 8.3 13.7 13.2 ...

summary(seeboden_clean)

##      stn                 time               tre200nn     
##  Length:11889       Min.   :1991-06-14   Min.   :-20.10  
##  Class :character   1st Qu.:1999-08-03   1st Qu.:  0.20  
##  Mode  :character   Median :2007-09-22   Median :  5.50  
##                     Mean   :2007-09-22   Mean   :  5.25  
##                     3rd Qu.:2015-11-11   3rd Qu.: 10.60  
##                     Max.   :2023-12-31   Max.   : 24.40  
##                                          NA's   :119

2.2.1 Visualizing

tre200nn = Lufttemperatur 2 m ¸ber Boden; Nachtminimum (18 UTC Vortag bis 6 UTC aktueller Tag)

#### EDA ####
# Histogram
hist_engelberg <- ggplot(engelberg_clean, aes(x = tre200nn)) +
  geom_histogram(bins = 30, fill = "blue", color = "black") +
  labs(title = "Distribution of Engelberg Temperature", x = "Temperature", y = "Count")

hist_seeboden <- ggplot(seeboden_clean, aes(x = tre200nn)) +
  geom_histogram(bins = 30, fill = "orange", color = "black") +
  labs(title = "Distribution of Seebodenalp Temperature", x = "Temperature", y = "Count")

# Boxplot
boxplot_engelberg <- ggplot(engelberg_clean, aes(y = tre200nn)) +
  geom_boxplot(fill = "blue", color = "darkred") +
  labs(title = "Boxplot of Engelberg Temperature", y = "Temperature")

boxplot_seeboden <- ggplot(seeboden_clean, aes(y = tre200nn)) +
  geom_boxplot(fill = "orange", color = "darkred") +
  labs(title = "Boxplot of Seebodenalp Temperature", y = "Temperature")

# Timeseries
timeseries_engelberg <- ggplot(engelberg_clean, aes(x = time, y = tre200nn)) +
  geom_line() +
  labs(title = "Time Series of Engelberg Temperature", x = "Time", y = "Temperature")

timeseries_seeboden <- ggplot(seeboden_clean, aes(x = time, y = tre200nn)) +
  geom_line() +
  labs(title = "Time Series of Seebodenalp Temperature", x = "Time", y = "Temperature")

2.2.1.1 Histogram

hist_engelberg

This histogram shows the temperature distribution in Engelberg, with most occurrences between 0 and 5 degrees Celsius. The data spans a broad range from below -20 to above 20 degrees Celsius, showing significant variability in temperatures. The distribution appears slightly left-skewed, suggesting that lower temperatures are more frequent than higher ones.

hist_seeboden

This histogram displays the temperature distribution in Seebodenalp (Rigi), with a symmetric shape centered around 5 to 10 degrees Celsius, indicating a normal distribution pattern. The temperatures range from below -20 to above 20 degrees Celsius, with the highest frequency of occurrences in the previously mentioned central range.

Compared to the Engelberg histogram, the Seebodenalp histogram shows a more normal distribution with temperatures clustered around a higher central range, whereas Engelberg’s distribution was slightly left-skewed with a mode closer to 0 to 5 degrees Celsius. Additionally, Seebodenalp’s temperature variability appears narrower with fewer extreme cold temperatures.

2.2.1.2 Boxplot

boxplot_engelberg

This boxplot of Engelberg temperature shows that the median temperature is slightly above 0 degrees Celsius, and the interquartile range (IQR) spans approximately from 5 to 15 degrees Celsius, indicating that the middle 50% of the temperatures are within this range. The whiskers extend from near -20 degrees Celsius to just above 20 degrees Celsius, suggesting that the majority of temperatures fall within this range, with some outliers on the lower end. The presence of several individual points below the lower whisker indicates outlier temperatures that are significantly colder than the rest.

boxplot_seeboden

The boxplot for Seebodenalp (Rigi) temperature indicates a median temperature close to 10 degrees Celsius with an interquartile range from approximately 0 to 20 degrees Celsius, suggesting that 50% of the temperatures lie within this warmer range. Similar to Engelberg, there are outliers shown below the bottom whisker, indicating occasional extreme cold temperatures.

Compared to Engelberg’s boxplot, the Seebodenalp boxplot demonstrates a higher median temperature and a narrower interquartile range, indicating less variability in temperature. Furthermore, while both locations show outliers, Engelberg appears to have a broader range of extreme cold temperatures.

2.2.1.3 Time series

timeseries_engelberg

What we see in this timeseries, is that there are missing data in the 1980s and it shows a seasonality. This is because the data still includes the summer months.

timeseries_seeboden

The data starts from the 1990s and also includes a seasonality aspect as the Engelberg timeseries does. This time series also presents higher temperature than the Engelberg one.

3 Time Series

library(tidyverse) # for data wrangling
library(lubridate) # date manipulation
library(TSstudio) # time series interactive viz
library(tseries) # for adf.test
library(astsa)
library(imputeTS)
library(forecast)
library(magrittr)

#### Clean Files ####

#homogenisierte Daten
Engelberg_homogenisiert <- read.csv("Daten/Messreihe_Engelberg.csv")
#temp_homo_ts <- ts(summary_data$Avg_Temperature, start = c(1864, 1), frequency = 12)

# Engelberg einlesen
engelberg <- read.csv("Daten/Engelberg_data.csv",sep=';', na.strings=c('-',''))

engelberg$time <- as.Date(as.character(engelberg$time), format = "%Y%m%d")
engelberg_clean <- engelberg %>% select('time',
                                        'tre200nn') %>% 
  rename('temp' = 'tre200nn') %>%
  filter(year(time) > 1989)

# Where are they missing
#Get dates for which temperatures are missing
missingCases <- which(is.na(engelberg_clean$temp)==TRUE)
u <- engelberg_clean$time[missingCases]

engelberg_clean <- engelberg_clean %>%
  filter(year(time) > 1989) %>%
  na_replace(fill=0)

ts_engelberg <- ts(data = engelberg_clean$temp,
               start = c(1990,01,01),
               frequency = 365)

3.1 Engelberg Data only

3.1.1 Trend Visualizing

# > Tägliche Daten von 1990-2023
engelberg_clean_trend <- lm(engelberg_clean$temp~engelberg_clean$time)

plot(ts_engelberg)
abline(engelberg_clean_trend, col = 'red')

3.1.2 Excursion: Monthly Data

We take only monthly data, which has been selected via a randomized procedure, to see how the Meteo and homogenous Data look like together.

# >> Monatliche Daten von 1990-2023 mit Trendlinie von Homogen.
temp_yr <- engelberg_clean %>%
  mutate(temp_raw=replace_na(temp,0)) %>%
  group_by(Year=year(time)) %>%
  filter(Year>=1990 & Year<=2023) %>%
  summarize(temp_yr=mean(temp_raw)) %>%
  ungroup()
temp_mn <- engelberg_clean %>%
  mutate(temp_raw=replace_na(temp,0)) %>%
  group_by(Year=year(time), Month=month(time)) %>%
  filter(Year>=1990 & Year<=2023) %>%
  summarize(temp_mn=mean(temp_raw), .groups='keep') %>%
  ungroup()

temp_mn_ts <- ts(temp_mn$temp_mn, start=c(temp_mn$Year[1],temp_mn$Month[1]), frequency=12)

plot(temp_mn_ts)

# Visualisierung des Trendes
fresh_snow_trend <- lm(temp_mn$temp_mn~temp_mn$Year)
engelberg_homog_trend <- lm(Engelberg_homogenisiert$Temperature~Engelberg_homogenisiert$Year)
plot(temp_mn_ts, xlab='Year', ylab='Average Temp.')
abline(fresh_snow_trend, col = 'red')
abline(engelberg_homog_trend, col = 'blue')

We will be exploring how the time series behaves for the whole year data.

# Yearly Data decomposing
ts_engelberg_dc <- decompose(ts_engelberg)
plot(ts_engelberg_dc)

# Stationary?
adf.test(ts_engelberg)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_engelberg
## Dickey-Fuller = -7.2169, Lag order = 22, p-value = 0.01
## alternative hypothesis: stationary

## it is stationary

acf(ts_engelberg)

acf(ts_engelberg_dc$random,na.action=na.pass)

# Plot the PACF
pacf(ts_engelberg_dc$random, na.action = na.pass)

# Plotting
PAutoCorrelation <- pacf(ts_engelberg_dc$random, na.action = na.pass, plot=FALSE)
plot(PAutoCorrelation, main = "Whole Year PACF")

#Arima
# # Fitting to PACF 
window <- list(start=1990,end=2023)

temp_comp_random_y <- data.frame(Date=time(ts_engelberg_dc$random), Random=ts_engelberg_dc$random)
temp_comp_random_y %<>% filter(Date>=window$start&Date<window$end)
temp_comp_random_y_ts <- ts(temp_comp_random_y$Random)
arima(temp_comp_random_y_ts, c(1,0,0))

## 
## Call:
## arima(x = temp_comp_random_y_ts, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.6528    -0.0089
## s.e.  0.0070     0.0740
## 
## sigma^2 estimated as 7.745:  log likelihood = -28625.83,  aic = 57257.65

The result of this Arima-Model is not a good simulation. Therefore, to focus on our core question: How many days are the Ski resorts able to use the artificial snow, we will be using only the winter months.

3.1.3 Winter

This further investigation only takes into account how many days have there been in a month where the temperature dropped below 0°.

engelberg_filtered <- engelberg_clean %>%
  # Assuming 'time' is already a Date object. If not, convert it first with as.Date()
  filter(year(time) > 1989) %>%
  filter(month(time) %in% c(12, 1, 2, 3)) 

# below 0 per month
monthly_below_zero <- engelberg_filtered %>%
  filter(temp < 0) %>%
  mutate(year = year(time),
         month = month(time)) %>%
  group_by(year, month) %>%
  summarise(days_below_zero = n(), .groups = "drop") %>% # Drop the grouping
  mutate(year_month = paste(year, month, sep="-")) %>% # Create year-month column
  select(year_month, days_below_zero) # Select only the columns you want

monthly_below_zero <- monthly_below_zero %>%
  mutate(year_month = as.Date(paste(year_month, "01", sep="-")))

ts_month <- ts(data = monthly_below_zero$days_below_zero,
                   start = c(1990,01),
                   frequency = 4)

autoplot(ts_month)

ts_month_dc <- decompose(ts_month)
plot(ts_month_dc)

adf.test(ts_month)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_month
## Dickey-Fuller = -4.4889, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

The decomposed time series of monthly days that measured below 0° is stationary.

acf(ts_month)

acf(ts_month_dc$random,na.action=na.pass)

pacf(ts_month_dc$random,na.action=na.pass)

# Arima Model 
month_auto <- auto.arima(ts_month, D=1, d=1)
month_auto

## Series: ts_month 
## ARIMA(3,1,0)(2,1,0)[4] 
## 
## Coefficients:
##           ar1      ar2      ar3     sar1     sar2
##       -0.6536  -0.4609  -0.1771  -0.8653  -0.4880
## s.e.   0.0884   0.1018   0.1044   0.0903   0.0797
## 
## sigma^2 = 33.01:  log likelihood = -404.83
## AIC=821.67   AICc=822.36   BIC=838.78

month_auto.1 <- auto.arima(ts_month)
month_auto.1

## Series: ts_month 
## ARIMA(0,0,2)(2,0,0)[4] with non-zero mean 
## 
## Coefficients:
##          ma1      ma2    sar1    sar2    mean
##       0.0949  -0.2114  0.2002  0.1514  22.120
## s.e.  0.0883   0.0943  0.0882  0.0980   0.591
## 
## sigma^2 = 26.93:  log likelihood = -405.4
## AIC=822.81   AICc=823.47   BIC=840.15

Forecast of our Arima Model

# Forecasting
f <- forecast(month_auto, level=c(95), h=5*4)


# Plotting the forecast
plot(f, main = "Forecast for the Next 5 Winter years", 
     xlab = "Time", ylab = "Forecast")

# Assuming 'f' is a forecast object that has a 'mean' component
y_lower <- min(-10, min(f$mean))
y_upper <- 60

plot(f, main = "Forecast for the Next 5 Winter years", 
     xlab = "Time", ylab = "Forecast",
     ylim = c(y_lower, y_upper))

Now we will be analyzing daily data of the winter months

ts_daily <- ts(data = engelberg_filtered$temp,
               start = c(1990,01),
               frequency = 365)
autoplot(ts_daily)

# Decompose
ts_daily_dc <- decompose(ts_daily)
plot(ts_daily_dc)

adf.test(ts_daily)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_daily
## Dickey-Fuller = -12.447, Lag order = 15, p-value = 0.01
## alternative hypothesis: stationary

acf(ts_daily)

acf(ts_daily_dc$random,na.action=na.pass)

PAutoCorrelation_m <- pacf(ts_daily_dc$random,na.action=na.pass, plot=FALSE)
plot(PAutoCorrelation_m, main = "Winter month PACF")

# Arima Model
# # Fitting to PACF 
window <- list(start=1990,end=2023)

temp_comp_random <- data.frame(Date=time(ts_daily_dc$random), Random=ts_daily_dc$random)
temp_comp_random %<>% filter(Date>=window$start&Date<window$end)
temp_comp_random_ts <- ts(temp_comp_random$Random)

arima(temp_comp_random_ts, c(1,0,0))

## 
## Call:
## arima(x = temp_comp_random_ts, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.6737     0.0204
## s.e.  0.0122     0.1619
## 
## sigma^2 estimated as 10.25:  log likelihood = -9473.05,  aic = 18952.09

Interpretation: The ARIMA model output shown indicates that the best-fitting model for the daily temperature data for a winter location during the winter months is an ARIMA(1,0,0), also known as an AR(1) model. The coefficient for `ar1` is 0.6737, suggesting a positive autocorrelation where a higher temperature on one day is likely to be followed by a higher temperature the next day. The intercept of 0.0204 suggests a very small upward trend in the temperature data, although its standard error is quite large (0.1619) relative to the coefficient value, indicating that this trend is not statistically significant.

3.1.4 One Year Forecast

In the following chapter we will be exploring the winter months of 2000 to maybe see how the forecast will be. 2000 has been chosen with a heuristic process because it is not the warmest winter nor the coldest.

one_engelberg <- engelberg_filtered %>% 
  filter(year(time) == 2000)
         # > 1999 &
         #   year(time) < 2003)

two_engelberg <- engelberg_filtered %>% 
  filter(year(time) == 2001)

ts_one <- ts(data = one_engelberg$temp, frequency = 30)
autoplot(ts_one)

ts_two <- ts(data = two_engelberg$temp, frequency = 30)

ts_one_dc <- decompose(ts_one)
plot(ts_one_dc)

adf.test(ts_one)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_one
## Dickey-Fuller = -2.9608, Lag order = 4, p-value = 0.1775
## alternative hypothesis: stationary

ts_one_diff <- na.omit(diff(ts_one))
adf.test(ts_one_diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_one_diff
## Dickey-Fuller = -6.0305, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

ts_one_diff_dc <- decompose(ts_one_diff)
plot(ts_one_diff_dc)

acf(ts_one_diff)

acf(ts_one_diff_dc$random,na.action=na.pass)

PAutoCorrelation_m <- pacf(ts_one_diff_dc$random,
                           na.action=na.pass, plot=FALSE)
plot(PAutoCorrelation_m, main = "1 winter PACF 2000")

window <- list(start=2000,end=2002)

temp_comp_random_one <- data.frame(Date=time(ts_one_dc$random), Random=ts_one_dc$random)
#temp_comp_random_one %<>% filter(Date>=window$start&Date<window$end)
temp_comp_random_one_ts <- ts(temp_comp_random_one$Random)

winter_model_dc <- arima(temp_comp_random_ts, c(1,0,0))

summary(winter_model_dc)

## 
## Call:
## arima(x = temp_comp_random_ts, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.6737     0.0204
## s.e.  0.0122     0.1619
## 
## sigma^2 estimated as 10.25:  log likelihood = -9473.05,  aic = 18952.09
## 
## Training set error measures:
##                         ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.0008668858 3.201432 2.455616 -37.91163 360.3382 0.9182744
##                     ACF1
## Training set -0.02672547

This model is using the decomposed data. The ARIMA(1,0,0) model summary for daily temperatures during winter months shows a significant positive autocorrelation in temperature data, as indicated by the ar1 coefficient of 0.6737. Error metrics such as the Root Mean Squared Error (RMSE) of 3.201432 and Mean Absolute Error (MAE) of 2.455616 suggest the model’s predictions are reasonably close to the actual temperatures, although the Mean Percentage Error (MPE) of -37.91163 indicates a systematic underestimation by the model.

winter_model <- arima(ts_one, c(2,0,0))
winter_model

## 
## Call:
## arima(x = ts_one, order = c(2, 0, 0))
## 
## Coefficients:
##          ar1     ar2  intercept
##       0.6454  0.0351    -2.7220
## s.e.  0.0906  0.0910     0.9538
## 
## sigma^2 estimated as 11.7:  log likelihood = -323.46,  aic = 654.91

The ARIMA(2,0,0) model for the daily temperature of a winter location indicates a primary autoregressive component at lag 1 (ar1 coefficient of 0.6454) and a much smaller effect at lag 2 (ar2 coefficient of 0.0351), suggesting that past temperatures have a significant influence on future temperatures, with the most recent day having the strongest effect. The negative intercept of -2.7220 may imply a baseline adjustment from the mean temperature, but its practical significance should be cautiously interpreted due to its relatively large standard error of 0.9538 compared to the coefficient value.

f <- forecast(winter_model, level=c(95), h=length(ts_two))
plot(f)

# Now add the actual data points.
#lines(ts_two, col="red")
forecast_start <- length(f$model$x)
# lines(seq(forecast_start, forecast_start + length(ts_two) - 1), 
# ts_two, col="red")

The plot displays the forecast from an ARIMA(2,0,0) model with a non-zero mean for future values of a time series, likely representing daily temperatures. The forecast shows a leveling off of the temperature values, with the shaded area representing a 95% confidence interval indicating where future observations are likely to fall, with the uncertainty increasing as the forecast extends further into the future.

3.2 Seebodenalp (Rigi)

kuessnacht <- read.csv("Daten/Seebodenalp_data.csv",sep=';', na.strings=c('-',''))

kuessnacht$time <- as.Date(as.character(kuessnacht$time), format = "%Y%m%d")
kuessnacht_clean <- kuessnacht %>% select('time',
                                        'tre200nn') %>% 
  rename('temp' = 'tre200nn') %>%
  filter(year(time) > 1989)
kuessnacht_clean <- na.omit(kuessnacht_clean[!is.na(kuessnacht_clean$temp), ])

3.2.1 Winter

kuessnacht_filtered <- kuessnacht_clean %>%
  # Assuming 'time' is already a Date object. If not, convert it first with as.Date()
  filter(year(time) > 1989) %>%
  filter(month(time) %in% c(12, 1, 2, 3)) 

Now we will be focusing on monthly days below 0°

monthly_below_zero <- kuessnacht_filtered %>%
  filter(temp < 0) %>%
  mutate(year = year(time),
         month = month(time)) %>%
  group_by(year, month) %>%
  summarise(days_below_zero = n(), .groups = "drop") %>% # Drop the grouping
  mutate(year_month = paste(year, month, sep="-")) %>% # Create year-month column
  select(year_month, days_below_zero) # Select only the columns you want

monthly_below_zero_k <- monthly_below_zero %>%
  mutate(year_month = as.Date(paste(year_month, "01", sep="-")))

ts_month_k <- ts(data = monthly_below_zero_k$days_below_zero,
               start = c(1990,01),
               frequency = 4)
autoplot(ts_month_k)

# Decompose
ts_month_k_dc <- decompose(ts_month_k)
plot(ts_month_k_dc)

adf.test(ts_month_k)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_month_k
## Dickey-Fuller = -4.0566, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

acf(ts_month_k)

acf(ts_month_k_dc$random,na.action=na.pass)

pacf(ts_month_k_dc$random,na.action=na.pass)

month_auto_k <- auto.arima(ts_month_k, D=1, d=1)
month_auto_k

## Series: ts_month_k 
## ARIMA(0,1,1)(0,1,2)[4] 
## 
## Coefficients:
##           ma1     sma1    sma2
##       -0.8398  -1.2282  0.2582
## s.e.   0.0624   0.1099  0.0862
## 
## sigma^2 = 29.82:  log likelihood = -392.57
## AIC=793.14   AICc=793.47   BIC=804.42

The ARIMA model summary indicates that the best-fitting model for the time series `ts_month_k` is an ARIMA(0,1,1)(0,1,2)[4], which suggests a seasonal model with a non-seasonal differencing of 1 and a seasonal differencing of 1 at a seasonal period of 4. The coefficients for the moving average part (ma1) and the first two seasonal moving average parts (sma1 and sma2) are significant, with ma1 being negative, indicating a smoothing effect from the previous value, and sma1 being strongly negative, which may suggest a corrective effect from the past seasonal cycle. The positive sma2 indicates a smaller smoothing effect from two seasons ago. The relatively low AIC, AICC, and BIC values suggest a model that balances goodness of fit with complexity, although model validation against unseen data is needed to confirm its predictive power.

month_auto_k.1 <- auto.arima(ts_month_k)
month_auto_k.1

## Series: ts_month_k 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     mean
##       -0.5957  0.7870  17.0877
## s.e.   0.1479  0.1077   0.5711
## 
## sigma^2 = 34.39:  log likelihood = -409.78
## AIC=827.55   AICc=827.87   BIC=838.99

f <- forecast(month_auto_k, level=c(95), h=5*4)

# Assuming 'f' is a forecast object that has a 'mean' component
y_lower <- min(-10, min(f$mean))
y_upper <- 60

plot(f, main = "Forecast for the Next 5 Winter", 
     xlab = "Time", ylab = "Forecast",
     ylim = c(y_lower, y_upper))

We will not be doing any other Time Series Analysis as we did for Engelberg, because for our analysis it is only of importance the below monthly 0° days for the winter months.

3.3 Tourism Data

EB <- read.csv("Daten/Room_Occupancy_Egelberg.csv")
KS <- read.csv("Daten/Room_Occupancy_Kuessnacht SZ.csv")

# Create time series
freq_monthly <- 12
Occupancy_Engelberg<-
  ts(EB$Room.occupancy, start=c(year(min(EB$Date)),yday(min(EB$Date))), frequency=freq_monthly) %>%
  na_replace(fill=0)
Occupancy_Kuessnacht<-
  ts(KS$Room.occupancy, start=c(year(min(KS$Date)),yday(min(KS$Date))), frequency=freq_monthly) %>%
  na_replace(fill=0)

plot(Occupancy_Engelberg)

plot(Occupancy_Kuessnacht)

# Sationarity test and decomposition
adf.test(Occupancy_Engelberg,k=0)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Occupancy_Engelberg
## Dickey-Fuller = -5.5441, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary

adf.test(Occupancy_Engelberg,k=20)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Occupancy_Engelberg
## Dickey-Fuller = -2.8082, Lag order = 20, p-value = 0.2408
## alternative hypothesis: stationary

Occupency_Engelberg_comp=decompose(Occupancy_Engelberg)
plot(Occupency_Engelberg_comp)

Occupency_Kuessnacht_comp=decompose(Occupancy_Kuessnacht)
plot(Occupency_Kuessnacht_comp)

adf.test(Occupancy_Engelberg,k=3)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Occupancy_Engelberg
## Dickey-Fuller = -4.552, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

acf(Occupency_Engelberg_comp$random,na.action=na.pass)

pacf(Occupency_Engelberg_comp$random,na.action=na.pass)

pacf(Occupancy_Engelberg)

adf.test(Occupancy_Kuessnacht,k=3)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Occupancy_Kuessnacht
## Dickey-Fuller = -7.1425, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

acf(Occupency_Kuessnacht_comp$random,na.action=na.pass)

pacf(Occupency_Kuessnacht_comp$random,na.action=na.pass)

pacf(Occupancy_Kuessnacht)

plot(Occupency_Engelberg_comp$random)

plot(Occupency_Kuessnacht_comp$random)

3.3.1 Cross Correlation Engelberg

ccf(na.omit(Occupency_Engelberg_comp$random), 
    na.omit(Occupency_Kuessnacht_comp$random))

Remove all non-winter months

# Delete all months except winter
EB$Date <- as.Date(EB$Date)
df_filtered <- EB[format(EB$Date, "%m") %in% c("12", "01", "02", "03"), ]
KS$Date <- as.Date(KS$Date)
df2_filtered <- KS[format(KS$Date, "%m") %in% c("12", "01", "02", "03"), ]

#########################################################################################3
# Time series temperature Engelberg
df_temp_EN <- read.csv("Daten/Engelberg_monthly_below_zero.csv")

df_temp_EN <- slice(df_temp_EN, -(1:92))

df_filtered <- slice(df_filtered, -42)
freq_monthly <- 4
Occupancy_Engelberg_season<-
  ts(as.numeric(df_filtered$Room.occupancy), frequency=freq_monthly) %>%
  na_replace(fill=0)
temp_Engelberg <- 
  ts(as.numeric(df_temp_EN$days_below_zero), frequency=freq_monthly) %>%
  na_replace(fill=0)

plot(Occupancy_Engelberg_season)

plot(temp_Engelberg)

Occupancy_Engelberg_season_comp <- decompose(Occupancy_Engelberg_season)
temp_Engelberg_comp <- decompose(temp_Engelberg)

plot(Occupancy_Engelberg_season_comp)

plot(temp_Engelberg_comp)

ccf(temp_Engelberg_comp$random, Occupancy_Engelberg_season_comp$random, na.action = na.pass)

3.3.2 Cross Correlation Küssnacht

df_temp_KS <- read.csv("Daten/kuessnacht_monthly_below_zero.csv")

df_temp_KS <- head(df_temp_KS, -2)

df_temp_KS <- tail(df_temp_KS, 42)

freq_monthly <- 4
Occupancy_Kuessnacht_season<-
  ts(as.numeric(df2_filtered$Room.occupancy), frequency=freq_monthly) %>%
  na_replace(fill=0)
temp_Kuessnacht <- 
  ts(as.numeric(df_temp_KS$days_below_zero), frequency=freq_monthly) %>%
  na_replace(fill=0)

plot(Occupancy_Kuessnacht_season)

plot(temp_Kuessnacht)

Occupancy_Kuessnacht_season_comp <- decompose(Occupancy_Kuessnacht_season)
temp_Kuessnacht_comp <- decompose(temp_Kuessnacht)
plot(Occupancy_Kuessnacht_season_comp)

plot(temp_Kuessnacht_comp)

ccf(temp_Kuessnacht_comp$random, Occupancy_Kuessnacht_season_comp$random, na.action = na.pass)

4 Conclusion

In conclusion, the analysis underscores the necessity for ski resorts to leverage advanced technologies and strategic diversification to adapt climate variability and ensure long-term viability. Using AI-driven weather prediction tools is essential in order to plan snowmkaing and maintenance, optimizing resort operations while complying with environmental regulations.

Additionally, ski resorts must continuously review and adjust operational practices to anticipate regulatory changes. By staying proactive, resorts can adapt to legal requirements, mitigate environmental impacts, and maintain operational efficiency.

Another important aspect is the diversification of resort attractions, with investments in summer sports infrastructure to ensure year-round tourism resilience. This approach mitigates the risk associated with seasonal dependency and contributes to economic stability. However, the ability to predict snow and weather conditions is limited due to a lack of extensive date from weather stations. Expanding the network of such stations across ski resorts could substantially improve this situation. By strategic placing, taking into account factors such as whether ski trails lie in sunny or shady areas, a deeper understanding of local climate variations could be gained. Such detailed climate data would significantly refine weather predictions and snowmaking plans, leading to more efficient and effective use of resources in resort operations.

Considering the research question about predicting the duration of effective snow-making using weather trends, our current predictive capabilities are limited. The scarcity of detailed weather data and the complex nature of ski trail environments make accurate forecasting challenging. Given additional time for investigation, we would focus on collecting more accurate meteorological data, taking into account the particular characteristics of ski trail locations such as altitude, exposure to sunlight, and patterns of shade. This precise approach to data gathering would enable us to develop more accurate predictive models and provide informed forecasts about the future of snowmaking and its impact on winter sports destinations.